***=====================CONTENT===========================
***Figure 1: Air Pollution and the Number of Births in China
    ***(a) Three-year Change in PM2.5
    ***(b) Yearly Mean PM2.5 

***Figure 2: Air Pollution and the Age-Specific Fertility Rate
    ***(a) Age 20-25, Three-Year PM2.5 Change  
    ***(b) Age 20-25, Yearly PM2.5 
	***(c) Age 26-30, Three-Year PM2.5 Change
    ***(d) Age 26-30, Yearly PM2.5
	
***Figure 3: Air Pollution and the Number of Ideal Children 
    ***(a) Three-year Change in PM2.5 
    ***(b) Yearly Mean PM2.5 
	
***Figure 4: Thermal Inversions and Air Pollution
    ***(a) The Strength of Thermal Inversion and PM2.5: 1998~2016
    ***(b) The Strength of Thermal Inversion and PM2.5: 1998~2016	
   
***Figure 5: Wind direction, distance, and coal consumption of large coal-fired plants 

***Figure A1: Thermal Inversions, Yearly PM2.5 and Three-year Changes in PM2.5

***Figure A2: Thermal Inversions, Elevation and PM2.5
   
***======================================================

*Figure 1
use "$workdata\census_han_reg",clear
collapse ever_give_birth pm25 dif_pm25_3,by(des_city)
*1a
preserve
xtile pm25_100 = dif_pm25_3, nq(100)
bysort pm25_100: egen pm25_mean_100=mean(dif_pm25_3)
bysort pm25_100: egen givebirth_100 = mean(ever_give_birth) 

set scheme s1color
twoway  scatter givebirth_100 pm25_mean_100, mcolor(green) msize(large) msymbol(Oh) mlwidth(medthick) ||   ///
        lfit    givebirth_100 pm25_mean_100,clwidth(thick) ||, ///
			   legend(off)   scale(1.0)    	 ///
ytitle("Num of births per reproductive women", size(med)) ///
xtitle("Three-year PM2.5 Change", size(med))
restore

*1b
xtile pm25_100 = pm25, nq(100)
bysort pm25_100: egen pm25_mean_100=mean(pm25)
bysort pm25_100: egen givebirth_100 = mean(ever_give_birth) 

set scheme s1color
twoway  scatter givebirth_100 pm25_mean_100, mcolor(green) msize(large) msymbol(Oh) mlwidth(medthick) ||   ///
        lfit    givebirth_100 pm25_mean_100,clwidth(thick) ||, ///
			   legend(off)   scale(1.0)    	 ///
ytitle("Num of births per reproductive women", size(med)) ///
xtitle("Annual Mean PM2.5", size(med)) ymtick(0.2(0.1)0.5) ylabel(0.2(0.1)0.5)

*Figure 2
use "$workdata\census_han_reg",clear
collapse ever_give_birth pm25 dif_pm25_3,by(des_city)
keep if age>20&age<=25

*2a
preserve
xtile pm25_100 = dif_pm25_3, nq(100)
bysort pm25_100: egen pm25_mean_100=mean(dif_pm25_3)
bysort pm25_100: egen givebirth_100 = mean(ever_give_birth) 

set scheme s1color
twoway  scatter givebirth_100 pm25_mean_100, mcolor(green) msize(large) msymbol(Oh) mlwidth(medthick) ||   ///
        lfit    givebirth_100 pm25_mean_100,clwidth(thick) ||, ///
			   legend(off)   scale(1.0)    	 ///
ytitle("Num of births per reproductive women", size(med)) ///
xtitle("Three-year PM2.5 Change", size(med))
restore

*2b
xtile pm25_100 = pm25, nq(100)
bysort pm25_100: egen pm25_mean_100=mean(pm25)
bysort pm25_100: egen givebirth_100 = mean(ever_give_birth) 

set scheme s1color
twoway  scatter givebirth_100 pm25_mean_100, mcolor(green) msize(large) msymbol(Oh) mlwidth(medthick) ||   ///
        lfit    givebirth_100 pm25_mean_100,clwidth(thick) ||, ///
			   legend(off)   scale(1.0)    	 ///
ytitle("Num of births per reproductive women", size(med)) ///
xtitle("Annual Mean PM2.5", size(med)) ymtick(0.2(0.1)0.5) ylabel(0.2(0.1)0.5)

*Figure 3
use "$workdata\clds_willing",clear
collapse dream_child_num pm25 dif_pm25_3,by(des_city)

*3a
preserve
xtile pm25_100 = dif_pm25_3, nq(100)
bysort pm25_100: egen pm25_mean_100=mean(dif_pm25_3)
bysort pm25_100: egen givebirth_100 = mean(dream_child_num) 

set scheme s1color
twoway  scatter givebirth_100 pm25_mean_100, mcolor(green) msize(large) msymbol(Oh) mlwidth(medthick) ||   ///
        lfit    givebirth_100 pm25_mean_100,clwidth(thick) ||, ///
			   legend(off)   scale(1.0)    	 ///
ytitle("Desired Num of Children", size(med)) ///
xtitle("Three-year PM2.5 Change", size(med))
restore

*3b
xtile pm25_100 = pm25, nq(100)
bysort pm25_100: egen pm25_mean_100=mean(pm25)
bysort pm25_100: egen givebirth_100 = mean(dream_child_num) 

set scheme s1color
twoway  scatter givebirth_100 pm25_mean_100, mcolor(green) msize(large) msymbol(Oh) mlwidth(medthick) ||   ///
        lfit    givebirth_100 pm25_mean_100,clwidth(thick) ||, ///
			   legend(off)   scale(1.0)    	 ///
ytitle("Desired Num of Children", size(med)) ///
xtitle("Annual Mean PM2.5", size(med))  ymtick(1(0.5)2.5) ylabel(1(0.5)2.5) 

*Figure 4
*(a)The Strength of Thermal Inversion and PM2.5: 1998~2016
use "$rawdata\Thermal", clear
g cityid = code/100
sort cityid year
merge cityid year using "$rawdata\pm25_1998_2017.dta"
keep if _merge==3
drop _merge
keep if year<2017

g thermalinversion_strength=abs(thermalstrength)
bysort year: egen mean_thermalinversion=mean(thermalinversion_strength)
bysort year: egen mean_pm25=mean(pm25)

gen group_num=(thermalinversion_strength>mean_thermalinversion)
gen group_thermalinversion=group_num*thermalinversion_strength
bysort year: egen total_thermalinversion=sum(thermalinversion_strength)
bysort year: egen upmean_thermalinversion=sum(group_thermalinversion)
gen lowmean_thermalinversion=total_thermalinversion-upmean_thermalinversion
twoway connected upmean_thermalinversion lowmean_thermalinversion year,scheme(s1color)

gen oneindex=1
bysort year: egen num_total=sum(oneindex)
bysort year: egen num_up_1=sum(group_num)
gen num_low_1=num_total-num_up_1

gen group1_pm25=group_num*pm25
bysort year: egen total_pm25=sum(pm25)
bysort year: egen up_1_pm25=sum(group1_pm25)
gen low_1_pm25=total_pm25-up_1_pm25
replace up_1_pm25=up_1_pm25/num_up_1
replace low_1_pm25=low_1_pm25/num_low_1

label var  mean_thermalinversion "The Strength of Thermal Inversion"
label var  up_1_pm25  "PM2.5(Strength of Thermal Inversion>  Mean Value)"
label var  low_1_pm25 "PM2.5(Strength of Thermal Inversion<= Mean Value)"
set scheme s1color
twoway ///
(line mean_thermalinversion year,  lp(longdash) yaxis(1))  ///
(connected up_1_pm25 low_1_pm25 year, yaxis(2) msymbol(circle triangle) mlwidth(thick thick))  ///
,title("The Strength of Thermal Inversion and PM2.5: 1998~2016") legend(cols(1))

*(b)The Strength of Thermal Inversion and PM2.5: 1998~2016
xtile quati=thermalstrength ,nq(100)
bysort quati: egen pm25_quati=mean(pm25)

twoway scatter pm25_quati quati, mcolor(green) msize(medlarge) msymbol(O) mlwidth(medthick)  ///
title("The Strength of Thermal Inversion and PM2.5: 1998~2016")  ///
xtitle("Quantile of Thermal Inversion(Strength)", size(med)) ///
ytitle("PM2.5",size(med))  ///
xmtick(0(10)100)  xlabel(0(10)100) ymtick(10(10)50) ylabel(10(10)50)

*Figure 5
use "$workdata\total_coal_2009",clear
merge 1:1 des_city using "$workdata\pm25_2009"
keep if _m == 3
drop _m
replace total_coal50_450 = total_coal50_450/1000
duplicates drop des_city,force
xtile quati=total_coal50_450,nq(100)
bysort quati: egen pm25_quati=mean(pm25)
bysort quati: egen coal_quati=mean(total_coal50_450)
set scheme s1color
twoway scatter pm25_quati coal_quati, mcolor(green) msize(large) msymbol(Oh) mlwidth(medthick) ||  ///
       lfit    pm25_quati coal_quati, clwidth(thick) ||,  ///
			   legend(off) scale(1.0)  ymtick(0(20)80) ylabel(0(20)80)  ///
xtitle("Wind IV", size(med)) ///
ytitle("PM2.5",size(med))  

*Figure A1
*a
use "$workdata/thermal_2009", clear
merge 1:1 des_city using "$workdata/pm25_2009"
keep if _m == 3
drop _m

xtile thermal_100 = thermalstrength2009, nq(100)
bysort thermal_100: egen pm25_mean_100=mean(pm25)
bysort thermal_100: egen thermal_mean_100 = mean(thermalstrength2009) 

set scheme s1color
twoway scatter pm25_mean_100 thermal_mean_100, mcolor(green) msize(large) msymbol(Oh) mlwidth(medthick) ||  ///
       lfit    pm25_mean_100 thermal_mean_100, clwidth(thick) ||,  ///
			   legend(off) scale(1.0)   ///
xtitle("The Strength of Thermal Inversion", size(med)) ///
ytitle("Yearly Mean PM2.5",size(med)) 

*b
use "$workdata/thermal_2009", clear
merge 1:1 des_city using "$workdata/pm25_2009"
keep if _m == 3
drop _m

xtile dif_thermal_100 = dif_thermalstrength_3, nq(100)
bysort dif_thermal_100: egen dif_pm25_mean_100=mean(dif_pm25_3)
bysort dif_thermal_100: egen dif_thermal_mean_100 = mean(dif_thermalstrength_3) 

set scheme s1color

twoway scatter dif_pm25_mean_100 dif_thermal_mean_100, mcolor(green) msize(large) msymbol(Oh) mlwidth(medthick) ||  ///
       lfit    dif_pm25_mean_100 dif_thermal_mean_100, clwidth(thick) ||,  ///
			   legend(off) scale(1.0)   ///
xtitle("Three-year Thermal Inversion Change", size(med)) ///
ytitle("Three-year PM2.5 Change",size(med)) 

*Figure A2
use "$workdata\elevation_wk",clear
egen mean_elevation = mean(alttud)
g high_city = alttud>mean_elevation
save "$workdata\elevation_wk_high",replace


use "$rawdata\Thermal", clear
g cityid = code/100
sort cityid year
merge cityid year using "$rawdata\pm25_1998_2017.dta"
keep if _merge==3
drop _merge

merge m:1 cityid using "$workdata\elevation_wk_high"
keep if _m==3
keep if year<2017

g thermalinversion_strength=abs(thermalstrength)
bysort year: egen mean_thermalinversion=mean(thermalinversion_strength)
bysort year: egen mean_pm25=mean(pm25)

gen group_num=(thermalinversion_strength>mean_thermalinversion)
gen group_thermalinversion=group_num*thermalinversion_strength
bysort year: egen total_thermalinversion=sum(thermalinversion_strength)
bysort year: egen upmean_thermalinversion=sum(group_thermalinversion)
gen lowmean_thermalinversion=total_thermalinversion-upmean_thermalinversion
twoway connected upmean_thermalinversion lowmean_thermalinversion year,scheme(s1color)

gen oneindex=1
bysort year: egen num_total=sum(oneindex)
bysort year: egen num_up_1=sum(group_num)
gen num_low_1=num_total-num_up_1

gen group1_pm25=group_num*pm25
bysort year: egen total_pm25=sum(pm25)
bysort year: egen up_1_pm25=sum(group1_pm25)
gen low_1_pm25=total_pm25-up_1_pm25
replace up_1_pm25=up_1_pm25/num_up_1
replace low_1_pm25=low_1_pm25/num_low_1

xtile quati=thermalstrength ,nq(100)
bysort quati: egen pm25_quati=mean(pm25)

set scheme s1color
twoway scatter pm25_quati quati if high_city==1, mcolor(green) msize(medlarge) msymbol(O) mlwidth(medthick) ///
title("Thermal Inversion and PM2.5(Elevation>Mean value)" )  ///
xtitle("Quantile of Thermal Inversion(Strength)", size(med)) ///
ytitle("PM2.5",size(med))  ///
xmtick(0(10)100)  xlabel(0(10)100) ymtick(10(10)50) ylabel(10(10)50) 


twoway scatter pm25_quati quati if high_city==0, mcolor(green) msize(medlarge) msymbol(O) mlwidth(medthick)  ///
title("Thermal Inversion and PM2.5(Elevation≤Mean value)" )  ///
xtitle("Quantile of Thermal Inversion(Strength)", size(med)) ///
ytitle("PM2.5",size(med))  ///
xmtick(0(10)100)  xlabel(0(10)100) ymtick(10(10)50) ylabel(10(10)50) 

